home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Turnbull China Bikeride
/
Turnbull China Bikeride - Disc 1.iso
/
ARGONET
/
PD
/
MATHS
/
RLAB
/
RLAB125.ZIP
/
!RLaB
/
toolbox
/
bessely
< prev
next >
Wrap
Text File
|
1994-09-23
|
7KB
|
249 lines
////-------------------------------------------------------------------//
// Syntax: bessely ( ALPHA, X )
// Description:
// Bessely is the bessel function of the second kind of order ALPHA.
// Either ALPHA or X may a vector in which case a vector result of
// the same dimension is returned. If both ALPHA and X are vectors
// a matrix result of dimension length(ALPHA) x length(X) is returned.
// If either ALPHA or X is a scalar, the other argument may be a matrix
// and a matrix result of the same dimension is returned.
//-------------------------------------------------------------------//
// See Numerical Recipes in C (second edition)
// Currently only integer order supported.
static (besselyn, bessely0, bessely1);
bessely = function(alph,x)
{
global (pi)
// Checking the class will automatically result in an error if
// the argument doesn't exist.
if(class(alph) != "num" || class(x) != "num") {
error("bessely: argument is non-numeric.");
}
if(min(size(alph)) > 1 && min(size(x)) > 1) {
// Both arguments are matrices.
error("bessely: only one argument may be a matrix.");
}
if(length(alph) == 1) {
// x is a matrix (or scalar), alph a scalar
if(alph == 0 || mod(alph,int(alph)) == 0) {
return besselyn(alph,x);
else
error("bessely: fractional orders not yet supported.");
// return besselyr(alph,x);
}
else if(length(x) == 1) {
// alph is a matrix (or vector), x a scalar
y = zeros(size(alph));
// Use integer order routines for integer alph
inte_ind = find(alph == 0 || mod(alph,int(alph)) == 0);
frac_ind = complement(inte_ind,1:alph.n);
if(inte_ind.n) {
for(index in inte_ind) {
y[index] = besselyn(alph[index],x);
}
}
if(frac_ind.n) {
error("bessely: fractional orders not yet supported.");
// for(index in frac_ind) {
// y[index] = besselyr(alph[index],x);
// }
}
return y;
else
// alph and x are both vectors
y = zeros(length(alph),length(x));
// Use integer order routines for integer alph
inte_ind = find(alph == 0 || mod(alph,int(alph)) == 0);
frac_ind = complement(inte_ind,1:alph.n);
if(inte_ind.n) {
for(index in inte_ind) {
y[index;] = besselyn(alph[index],x);
}
}
if(frac_ind.n) {
error("bessely: fractional orders not yet supported.");
// for(index in frac_ind) {
// y[index;] = besselyr(alph[index],x);
// }
}
return y;
}}
};
// In these functions x can be a vector (or matrix), but n
// must be a scalar.
// No argument checking is performed on static functions.
besselyn = function ( n, x ) {
local(abs_x,y,index,index_1,index_2,index_3,arg,p_1,p_2,p_3, ...
two_over_x,m,y_2,Norm,even,limit,lim_ind);
if(n == 0) {
return bessely0(x);
else if (n == 1) {
return bessely1(x);
else
// integer order greater than two
y = zeros(size(x));
index_1 = 1:x.n;
index_2 = find(x == 0);
index_1 = complement(intersection(index_1,index_2),index_1);
if(index_2.n) {
y[index_2] = -1./zeros(size(index_2));
}
if(index_1.n) {
arg = x[index_1];
p_1 = bessely0(arg);
p_2 = bessely1(arg);
two_over_x = 2./arg;
for(index in 1:(n-1)) {
p_3 = index * (two_over_x .* p_2) - p_1;
p_1 = p_2;
p_2 = p_3;
}
y[index_1] = p_3;
}
return y;
}}
};
bessely0 = function(x) {
local(index_1,index_2,index_3,y,c_1,c_2,arg_1,arg_2,p_1,p_2);
index_1 = find(x < 8);
index_2 = complement(index_1,1:x.n);
index_3 = find(x == 0.);
index_1 = complement(intersection(index_1,index_3),index_1);
// First evaluate using rational approx. for small arguments.
y = zeros(size(x));
if(index_3.n) {
y[index_3] = -1./zeros(size(index_3));
}
if(index_1.n) {
c_1 = [ -2957821389, 7062834065, -512359803.6, ...
10879881.29, -86327.92757, 228.4622733 ];
c_2 = [ 40076544269, 745249964.8, 7189466.438, ...
47447.26470, 226.1030244 ];
arg_1 = x[index_1];
arg_2 = arg_1 .* arg_1;
// Rational approximation to the "regular part"
p_1 = c_1[1] + arg_2 .* (c_1[2] + arg_2 .* (c_1[3] + arg_2 .* ...
(c_1[4] + arg_2 .* (c_1[5] + arg_2 .* c_1[6]))));
p_2 = c_2[1] + arg_2 .* (c_2[2] + arg_2 .* (c_2[3] + arg_2 .* ...
(c_2[4] + arg_2 .* (c_2[5] + arg_2))));
y[index_1] = p_1./p_2 + (2/pi)*besselj(0,arg_1).*log(arg_1);
}
// Now evaluate for large arguments using approximating form.
if(index_2.n) {
c_1 = [ -0.1098628627e-2, 0.2734510407e-4, -0.2073370639e-5, ...
0.2093887211e-6 ];
c_2 = [ -0.1562499995e-1, 0.1430488765e-3, -0.6911147651e-5, ...
0.7621095161e-6, -0.934945152e-7 ];
arg_1 = 8./x[index_2];
arg_2 = arg_1 .* arg_1;
p_1 = 1 + arg_2 .* (c_1[1] + arg_2 .* (c_1[2] + arg_2 .* (c_1[3] + ...
arg_2 .* c_1[4])));
p_2 = c_2[1] + arg_2 .* (c_2[2] + arg_2 .* (c_2[3] + arg_2 .* (c_2[4] + ...
arg_2 .* c_2[5])));
c_1 = x[index_2];
c_2 = c_1 - pi/4;
y[index_2] = sqrt(2./(pi*c_1)).*(p_1.*sin(c_2) + arg_1.*p_2.*cos(c_2));
}
return y;
};
bessely1 = function(x) {
local(index_1,index_2,index_3,y,c_1,c_2,arg_1,arg_2,p_1,p_2);
index_1 = find(x < 8);
index_2 = complement(index_1,1:x.n);
index_3 = find(x == 0.);
index_1 = complement(intersection(index_1,index_3),index_1);
// First evaluate using rational approx. for small arguments.
y = zeros(size(x));
if(index_3.n) {
y[index_3] = -1./zeros(size(index_3));
}
if(index_1.n) {
c_1 = [ -0.4900604943e13, 0.1275274390e13, -0.5153438139e11, ...
0.7349264551e9, -0.4237922726e7, 0.8511937935e4 ];
c_2 = [ 0.2499580570e14, 0.4244419664e12, 0.3733650367e10, ...
0.2245904002e8, 0.1020426050e6, 0.3549632885e3 ];
arg_1 = x[index_1];
arg_2 = arg_1 .* arg_1;
p_1 = arg_1.*(c_1[1] + arg_2 .* (c_1[2] + arg_2 .* (c_1[3] + arg_2 .* ...
(c_1[4] + arg_2 .* (c_1[5] + arg_2 .* c_1[6])))));
p_2 = c_2[1] + arg_2 .* (c_2[2] + arg_2 .* (c_2[3] + arg_2 .* ...
(c_2[4] + arg_2 .* (c_2[5] + arg_2 .* (c_2[6] +arg_2)))));
y[index_1] = p_1./p_2 + (2/pi)*(besselj(1,arg_1).*log(arg_1) - 1./arg_1);
}
// Now evaluate for large arguments using approximating form.
if(index_2.n) {
c_1 = [ 0.183105e-2, -0.3516396496e-4, 0.2457520174e-5, ...
-0.240337019e-6 ];
c_2 = [ 0.04687499995, -0.2002690873e-3, 0.8449199096e-5, ...
-0.88228987e-6, 0.105787412e-6 ];
arg_1 = 8./x[index_2];
arg_2 = arg_1 .* arg_1;
p_1 = 1 + arg_2 .* (c_1[1] + arg_2 .* (c_1[2] + arg_2 .* (c_1[3] + ...
arg_2 .* c_1[4])));
p_2 = c_2[1] + arg_2 .* (c_2[2] + arg_2 .* (c_2[3] + arg_2 .* (c_2[4] + ...
arg_2 .* c_2[5])));
c_1 = x[index_2];
c_2 = c_1 - 3*pi/4;
y[index_2] = sqrt(2./(pi*c_1)).*(p_1.*sin(c_2) + arg_1.*p_2.*cos(c_2));
}
return y;
};